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ABSTRACT 

The abundance and structure of dark matter subhalos has been analyzed extensively in recent 
studies of dark matter-only simulations, but comparatively little is known about the impact of 
baryonic physics on halo substructures. We here extend the SUBFIND algorithm for substruc- 
ture identification such that it can be reliably applied to dissipative hydrodynamical simula- 
tions that include star formation. This allows, in particular, the identification of galaxies as 
substructures in simulations of clusters of galaxies, and a determination of their content of 
gravitationally bound stars, dark matter, and hot and cold gas. Using a large set of cosmolog- 
ical cluster simulations, we present a detailed analysis of halo substructures in hydrodynam- 
ical simulations of galaxy clusters, focusing in particular on the influence both of radiative 
and non-radiative gas physics, and of non-standard physics such as thermal conduction and 
feedback by galactic outflows. We also examine the impact of numerical nuisance parame- 
ters such as artificial viscosity parameterizations. We find that diffuse hot gas is efficiently 
stripped from subhalos when they enter the highly pressurized cluster atmosphere. This has 
the effect of decreasing the subhalo mass function relative to a corresponding dark matter- 
only simulation. These effects are mitigated in radiative runs, where baryons condense in the 
central subhalo regions and form compact stellar cores. However, in all cases, only a very 
small fraction, of the order of one percent, of subhalos within the cluster virial radii preserve 
a gravitationally bound hot gaseous atmosphere. The fraction of mass contributed by gas in 
subhalos is found to increase with the cluster-centric distance. Interestingly, this trend extends 
well beyond the virial radii, thus showing that galaxies feel the environment of the pressurized 
cluster gas over fairly large distances. The compact stellar cores (i.e. galaxies) are generally 
more resistant against tidal disruption than pure dark matter subhalos. Still, the fraction of 
star-dominated substructures within our simulated clusters is only ^10 per cent. We expect 
that the finite resolution in our simulations makes the galaxies overly susceptible to tidal dis- 
ruption, hence the above fraction of star-dominated galaxies should represent a lower limit 
for the actual fraction of galaxies surviving the disruption of their host dark matter subhalo. 

Key words: hydrodynamics, method: numerical, galaxies: cluster: general, galaxies: evolu- 
tion, cosmology: theory 



1 INTRODUCTION 

The hierarchical cold dark matter (CDM) model has been estab- 
lished as the standard model of cosmic structure formation and 
its parameters are now tightly constrained by a variety of observa- 
tions (e.g. Komatsu et al. 2008; Springel et al. 2006 and references 
therein). Within the CDM scenario, galaxy clusters are the largest 
gravitationally bound objects in the universe, and they exhibit rich 
information about the process of structure formation. As a result, 
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they attract great interest both from observational and theoretical 
points of view. Due to their recent formation, clusters of galaxies 
also represent the ideal places where the hierarchical assembly of 
structures is 'caught in the act'. Indeed, the complexity of their in- 
ternal structure reflects the infall of hundreds, if not thousands, of 
smaller objects and their subsequent destruction or survival within 
the cluster potential. The dynamical processes determining the fate 
of the accreting structures also provide the link between the internal 
dynamics of clusters and the properties of their galaxy population. 

Modem cosmological simulations provide an ideal tool to 
study the non-linear dynamics of these processes in a cosmolog- 
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ical environment. Thanks to the high resolution reached by dark 
matter (DM) only simulations, a consistent picture of the popu- 
lation of subhalos within galaxy clusters has now emerged (e.g., 
Moore et al. 1999; Ghigna et al. 2000; Springel et al. 2001; Stoehr 
et al. 2002, 2003; Diemand et al. 2004; Gao et al. 2004; De Lucia 
et al. 2004; Kravtsov et al. 2004). Although these simulations give 
a highly detailed description of the evolution of DM substructures, 
they can not provide information on how gas-dynamical processes 
affect the properties and the dynamics of substructures. Yet such an 
understanding is needed to make a more direct link of the proper- 
ties of dark matter substructures with those of the galaxies that orbit 
in clusters. For example, we would like to know the timescale over 
which subhalos can retain their diffuse gas once they infall into a 
cluster, as this also determines the time interval over which gas can 
continue to cool and to refuel ongoing star formation. Being highly 
pressurized, the intra-cluster medivmi (ICM) is quite efficient in 
stripping gas from merging structures (Gunn & Gott 1972; see also 
McCarthy et al. 2008 and references therein), thereby altering both 
the mass of the subhalos and their survival time. The stripping also 
governs the timescale over which galaxy morphology and colours 
are modified inside galaxy clusters (Abadi et al. 1999; Kenney et al. 
2004). Finally, although baryons are expected to only play a sub- 
dominant role for the dynamics of substructures due to their limited 
mass fraction, they may nevertheless alter the structure of subhalos 
in important ways, affecting their survival, mass loss rate, abun- 
dance and radial distribution. Quantifying these differences relative 
to dark matter-only simulations is an important challenge for the 
present generation of cosmological hydrodynamical simulations. 

So far there have only been a limited nvunber of investiga- 
tions of the properties of sub-halos within hydrodynamical simula- 
tions. Based on cluster simulationsw with the adaptive mesh refine- 
ment code ART, Nagai & Kravtsov (2005) have demonstrated that 
the stellar mass of infalling galaxies is approximately conserved 
and the survival of the galaxies is slightly enhanced. Similar stud- 
ies using the Tree-SPH code GASOLINE have been performed by 
Maccio et al. (2006), where twice as many sub-halos within the 
central part of a Galaxy sized halo were found when compared with 
pure dark matter simulations. Finally Weinberg et al. (2008) used a 
TreeSPH code to investigate the halo occupation and sub-halo cor- 
relation functions for one group-sized halo, also pointing out the 
enhanced survival of sub-structures in the hydrodynamical runs. 

Here we present a systematic analysis of a large sample of 
galaxy clusters, with an emphasis on the role of non-standard 
physics and the details of feedback prescriptions on the survival of 
galaxies. Additionally, we also present an analysis of the evolution 
of the diffuse gas within the galaxies. In general our findings are in 
broad agreement with those in previous studies, but we stress that 
our star dominant systems are still not as numerous as expected 
from semi-analytic modelling to reproduce the cluster luminosity 
functions (e.g., De Lucia & Blaizot (2007); see also Baugh (2006) 
for a review on semi-analytic models of galaxy formation). 

Thanks to improvements in nimierical codes and higher com- 
puter performance, cosmological hydrodynamical simulations of 
galaxy clusters can now model a number of the most important 
physical processes in galaxy formation, such as radiative cooling, 
star formation, and feedback in energy and metals, while at the 
same time reaching sufficient resolution to treat these processes in 
a numerically robust and physically meaningful way (e.g., Borgani 
et al. 2008, for a recent review). This allows more realistic struc- 
ture formation calculations of the coupled system of dark matter 
and baryonic gas than possible with dark matter only simulations, 
far into the highly non-Unear regime. For instance, radiative cooUng 



has the effect of letting a significant fraction of the baryons cool at 
the halo centres, which alters the shape and concentration of halos 
through the mechanism of adiabatic contraction (e.g., Gnedin et al. 
2004). At the same time, the condensation of cooled baryons and 
their subsequent conversion into coUisionlcss stars is expected to 
largely protect them from stripping, reducing the rate of disruption 
of substructures and increasing their survival times. Indeed, assum- 
ing that galaxies at least survive for a while the disruption of their 
host DM halos has been shown to be crucial for semi-analytical 
models of galaxy formation to provide a successful description of 
the observed galaxy population in clusters (e.g., Gao et al. 2004). 

This discussion highlights the importance of understanding 
the properties and the evolution of cluster substructures as pre- 
dicted by modern cosmological hydrodynamical simulations. Ev- 
idently, an important technical prerequisite for such an analysis is 
the availability of suitable nimierical algorithms to reliably identify 
substructure in dissipative simulations. To this end we introduce a 
modified version of the SUBFIND algorithm, and test its robustness. 
We then apply it in a detailed analysis of the properties of substruc- 
tures in a large ensemble of galaxy clusters, which have been simu- 
lated both at different resolutions in their DM-only version, and in 
several further versions where different physical processes where 
included, such as non-radiative and radiative hydrodynamics, star 
formation, energy feedback, gas viscosity and thermal conduction. 
The aim of this analysis is to quantify the numerical robustness of 
the measured properties of substructures, and to identify the quan- 
titative influence of these physical processes on substructure statis- 
tics in comparison with dark matter only models. 

The paper is organized as follows. We describe in Section 2 the 
sample of galaxy clusters and the simulations performed. Section 3 
provides a description of the algorithm used to identify gravita- 
tionally bound substructures. This algorithm is a modified version 
of SUBFIND (Springel et al. 2001), which we suitably changed 
to allow extraction of bound structures from N-Body/SPH simula- 
tions also in the presence of gas and star particles. In Section 4, 
we present our results about the properties of the substructures in 
both DM and hydrodynamical simulations. We summarize our re- 
sults and outline our main conclusion in Section 5. An Appendix 
is devoted to the presentation of tests of resolution and numerical 
stabiUty of our analysis. 



2 SIMULATIONS 

In this section we describe our set of simulated clusters, the differ- 
ent physical processes we considered, and the numerical choices 
made for the different runs. 



2.1 Simulated physics 

The simulations were carried out with the TreePM/SPH code 
GADGET-2 (Springel et al. 2001; Springel 2005), which makes use 
of the entropy-conserving formulation of SPH (Springel & Hern- 
quist 2002). Some of our non-radiative simulations were carried 
out with different implementations of artificial viscosity compared 
with the standard one in GADGET-2, however, which allows us to 
investigate the effect of numerical viscosity on the stripping of gas 
from substructures within galaxy clusters. In Dolag et al. (2005) 
it was shown that the amount of turbulence detectable within the 
cluster atmosphere is quite sensitive to the numerical treatment of 
the artificial viscosity, so one expects that the stripping of the gas 
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Figure 1. Visualisation of the gas temperature in the non-radiative {ovisc) simulations of the ten clusters discussed in the text, car- 
ried out with the ray-tracing software SPLOTCH. An animation flying through one of the simulated clusters can be downloaded from 
http : / / www . mpa-garching . mpg . de/ galform/data_vis/ index . shtml#movie 9. 
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Name Included Physics Reference 


Used in 


dark matter only 


dm gravity - 
{g676, g914,gl542, g3344,g6212, gl, g8,g51, g72 and g696) 


Puchwein et al. (2005); Giocoli et al. (2008); Meneghetti et al. 
(2007, 2008) 


non radiative runs 


ovisc usual parametrisation of artificial Dolag et al. (2005) 
viscosity 

(g676, g914,gl542, g3344,g6212, gl, g8,g51, g72 and g696) 


Puchwein et al. (2005); Ettori et al. (2006); Dolag et al. (2006); 
Bonaldi et al. (2007); Cora et al. (2008); Puchwein & Bartel- 
mann (2006, 2007); Meneghetti et al. (2007) 


svisc artificial viscosity based on signal- Dolag et al. (2005) 
velocity 

(g676, g914,gl542, g3344,g6212, gl, g8,g51 and g72) 


Ivisc time varying, low artificial viscosity Dolag et al. (2005) 

scheme 

(g676, g914,gl542, g3344,g6212, gl, g8,g51 and g72) 


Ettori et al. (2006); Vazza et al. (2006) 


cooling, star-formation and thermal feedback 


csfnw no winds - 
(g676 and g51) 


Ettori et al. (2006); Borgani et al. (2006) 


csf weak winds (Vw = 340km/h) - 
(g676, g914,gl542, g3344,g6212, gl, g8,g51, g72 and g696) 


Puchwein et al. (2005); Rasia et al. (2005); Ettori et al. (2006); 
Bonaldi et al. (2007); Borgani et al. (2006); Ameglio et al. 
(2006, 2007) 


csfsw strong winds (vw = 480km/h) 
(g676 and g52) 


Ettori et al. (2006); Borgani et al. (2006) 


csfc weak winds {Vu, = 340kin/h) and Dolag et al. (2004) 

thermal conduction (k = 0.3) 
(g676, g914,gl542, g3344,g6212, gl, g8,g51, g72 and g696) 


Ettori et al. (2006); Rasia et al. (2006); Bonaldi et al. (2007) 



Table 1. Symbolic names of the different runs analysed here, together with a short description of the physical processes included, and references to studies 
that used the simulations previously or that give specific details for the physical models used. For each physical model, we also specify for which cluster the 
simulations have been carried out. 



from substructures might also be affected by these numerical de- 
tails. In the following, we will label simulations that use the origi- 
nal parametrisation of the artificial viscosity by Monaghan & Gin- 
gold (1983); Balsara (1995) as ovisc. This was found to be the 
scheme with the highest numerical viscosity in the study of Dolag 
et al. (2005). An alternative formulation with slightly less numeri- 
cal viscosity is based on the signal velocity approach of Monaghan 
(1997), and is labelled as svisc. Finally, we considered a modi- 
fied artificial viscosity scheme as originally suggested by Morris 
& Monaghan (1997), labelled as Ivisc. In this scheme, every par- 
ticle evolves its own time-dependent viscosity parameter, with the 
goal to only have high viscosity in regions where it is really needed. 
In this scheme, shocks are generally as well captured as in the stan- 
dard approach, but regions away from the shocks experience less ar- 
tificial viscosity, such that an inviscid ideal gas is represented more 
faithfully. As a result, turbulence driven by fluid instabilities can be 
better resolved, and simulated galaxy clusters are able to build up 
a higher level of turbulence generated along the shear flows arising 
in the cosmological structure formation process (see Dolag et al. 
2005). 

GADGET-2 also includes, if enabled, radiative cooling, heating 
by a uniform redshift-dependent UV backgroimd (Haardt & Madau 
1996), and a treatment of star formation and feedback processes. 



Simulations that account for radiative cooling and star formation 
allow us to study the effect of a compact stellar core at the centre of 
substructures on subhalo survival and disruption. The prescription 
of star formation we use is based on a sub-resolution model to 
account for the multi-phase structure of the interstellar medium 
(ISM), where the cold phase of the ISM is the reservoir of star 
formation (Springel & Hernquist 2003). Supernovae heat the hot 
phase of the ISM and provide energy for evaporating some of the 
cold clouds, thereby leading to self-regulation of the star formation 
and an effective equation of state for describing its dynamics. 

As a phenomenological extension of this feedback scheme, 
Springel & Hernquist (2003) also included a simple model for 
galactic winds, whose velocity, v-w, scales with the fraction r? of 
the SN type-n feedback energy that contributes to the winds. The 
total energy provided by SN type-ll is computed by assuming that 
they are due to exploding massive stars with mass > 8 M0 from a 
Salpeter (1955) initial mass function (IMF), with each SN releasing 
10^^ ergs of energy. In our simulations with winds, we will assume 
?7 = 0.5 and 1, yielding Vw ^ 340 (csf) and 480 km s"'^ (csfsw), 
respectively, while we will also explore the effect of switching off 
galactic winds (csfiiw). We refer to Table 1 for a schematic descrip- 
tion and overview of the physical and numerical schemes included 
in our simulations. 
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dm ovisc 



Cluster 






-R200 


M200 


-Rboo 




-Rvir 


iWvir 


.R200 


M200 


R500 


Mboo 


g676.a 


1.43 


1.60 


1.07 


1.39 


0.71 


1.06 


1.40 


1.53 


1.06 


1.33 


0.71 


1.03 


g914.a 


1.50 


1.84 


1.09 


1.44 


0.71 


1.06 


1.46 


1.73 


1.09 


1.43 


0.71 


1.06 


gl342.a 


1.40 


1 .53 


1 .043 


1.29 


0.69 


0.93 


1.40 




1.043 


1.30 


0.69 


0.94 


g3344.a 


1.44 


1.64 


1.07 


1.40 


0.73 


1.10 


1.43 


1.63 


1.07 


1.39 


0.73 


1.07 


g6212.a 


1.43 


1.61 


1.06 


1.33 


0.70 


1.00 


1.43 


1.61 


1.06 


1.31 


0.70 


1.00 


g51.a 


3.26 


19.21 


2.39 


15.33 


1.56 


10.66 


3.26 


19.03 


2.37 


15.30 


1.57 


11.16 


g72.a 
g72.b 


3.29 
1.71 


19.63 
2.81 


2.40 
1.21 


15.59 
2.03 


1.57 
0.71 


10.93 
1.04 


3.29 
1.71 


19.63 
2.79 


2.37 
1.26 


15.29 
2.24 


1.59 
0.74 


11.20 
1.19 



gl.a 


3.40 


21.86 


2.54 


18.86 


1.73 


14.61 


3.37 


21.16 


2.50 


17.80 


1.69 


13.59 


gl.b 


2.31 


6.83 


1.64 


5.07 


1.04 


3.26 


2.26 


6.43 


1.67 


5.36 


1.06 


3.43 


gl.c 


1.67 


2.61 


1.23 


2.09 


0.79 


1.40 


1.69 


2.66 


1.23 


2.09 


0.79 


1.41 


gl.d 


1.60 


2.29 


1.10 


1.53 


0.63 


0.71 


1.57 


2.14 


1.07 


1.40 


0.61 


0.66 


gl.e 


1.27 


1.14 


0.94 


0.94 


0.59 


0.56 


1.26 


1.10 


0.93 


0.93 


0.63 


0.70 


gl-f 


1.14 


0.81 


0.77 


0.53 


0.49 


0.34 


1.11 


0.77 


0.77 


0.53 


0.50 


0.34 


gS.a 


3.90 


32.70 


2.86 


26.63 


1.90 


19.64 


3.94 


33.87 


2.89 


27.56 


1.93 


20.60 


gS.b 


1.50 


1.86 


1.09 


1.44 


0.69 


0.93 


1.54 


2.01 


1.14 


1.70 


0.74 


1.20 


gS.c 


1.36 


1.37 


0.93 


0.93 


0.60 


0.60 


1.43 


1.61 


1.03 


1.24 


0.61 


0.66 


gS.d 


1.34 


1.34 


0.89 


0.80 


0.59 


0.56 


1.40 


1.51 


1.00 


1.14 


0.66 


0.81 


gS.e 


1 .30 


1.23 


0.96 


1.00 


0.60 


0.61 


1.33 


1.30 


0.96 


1.00 


0.63 


0.69 


gS.f 


1.14 


0.83 


0.83 


0.66 


0.46 


0.29 


1.20 


0.94 


0.87 


0.76 


0.54 


0.47 


gs.g 


1.11 


0.76 


0.77 


0.53 


0.47 


0.31 


1.06 


0.66 


0.76 


0.49 


0.50 


0.36 



g696.a 


3.26 


19.07 


2.40 


15.79 


1.56 


g696.b 


3.13 


17.04 


2.23 


12.64 


1.43 


g696.c 


2.79 


11.97 


2.07 


10.04 


1.27 


g696.d 


2.67 


10.57 


1.83 


7.04 


1.14 



10.81 


3.24 


18.94 


2.39 


15.51 


1.57 


11.00 


8.23 


3.19 


17.87 


2.29 


13.51 


1.47 


9.10 


5.81 


2.77 


11.87 


2.07 


10.13 


1.29 


6.06 


4.30 


2.59 


9.51 


1.81 


6.80 


1.14 


4.29 



Table 2. General properties of the simulated clusters at z = for the cases of DM-only runs (left part of the table) and for the ovisc version of the non-radiative 
runs (right part of the table). Column 1 : cluster names; cols. 2 (7) and 3 (8): virial radii (in units of Mpc) and virial masses (in units of 10^** M©); cols. 5 (10) 
and 4 (9): masses contained within the radius i?2oo. encompassing an average density of 200/Ocrit and values of R200 \ cols. 7 (12) and 6 (11): the same as 
before, but referring to a mean density of 500 Pcrit- 



For some of our cluster simulations we also included the effect 
of heat conduction (csfc), based on the implementation described 
by Jubelgas et al. (2004). In the simulations presented here, we as- 
sume an isotropic effective conductivity parameterized as a fixed 
fraction of 1/3 the Spitzer rate. We also account for saturation, 
which can become relevant in low-density gas and at the interface 
between cold and hot media. We refer to Dolag et al. (2004) for 
more details on the effect of thermal conduction on the thermody- 
namical properties of the intra-cluster medium (ICM). 



2.2 The set of simulated clusters 

The clusters analyzed in this study are extracted from 10 re- 
simulations of Lagrangian regions selected from a cosmological 
lower resolution DM-only simulation (Yoshida et al. 2001). This 
parent simulation has a box-size of 479 Mpc, and assumed 
a flat ACDM cosmology with Q,m = 0.3 for the matter den- 
sity parameter, Hq — 70 km s^^Mpc^^ for the Hubble constant, 
/bar = 0.13 for the baryon fraction and as = 0.9 for the normal- 
isation of the power spectrum. As such, the values of Q.m. and as 
in this model are somewhat higher than the best-fitting values ob- 
tained from the analysis of the 5-year WMAP data (Dunkley et al. 
2008), but are still consistent with the current cosmological con- 
straints. 



Five of the regions have been selected to surround low-mass 
clusters (mdustor ~ 10^* Mg), while four other regions surround 
massive (mciuster^ lO^^M©) clusters. Besides the central mas- 
sive cluster, three of these four regions also contain 12 additional 
smaller clusters, all having virial masses above 10^'^ Mq. The tenth 
simulated region surrounds a filamentary structure which includes 
four massive (mciuster > lO^''Af0) clusters (see Dolag et al. 
2006). Ray-tracing images, obtained with the SPLOTCH package 
(Dolag et al. 2008), of the gas temperature for the simulated clus- 
ter regions are shown in Figure 1. Besides giving an impression 
of the complexity and the dynamics of the ICM within the sim- 
ulated galaxy clusters, this figure also highlights a couple of in- 
teresting key features of hydrodynamical cluster simulations. The 
white colours in general correspond to high density, whereas the red 
colour marks the hot, shock-heated atmosphere of the cluster. The 
transition to the un-shocked material (mainly indicated by the blue 
colour) is clearly visible and indicates the location of the accretion 
shocks, which for clusters at redshift zero are typically at much 
larger distances than the virial radius (see Table 2). Also clearly 
visible are denser and colder filaments that penetrate the hot cluster 
atmosphere. Many of the gaseous halos of the substructures show 
comet like features, a tell tale sign of gas being stripped. Sometimes 
this happens already at relatively large distances from the cluster 
center, demonstrating a change in the enviroment that becomes im- 
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Figure 2. Visualisation of all the identified galaxies within the virial radius of the g51, g72, gl and gS clusters (from top left to bottom right), simulated 
with cooling and star-formation and galactic winds (csf). About 1000 galaxies are identified within the virial radius of each of these massive clusters. Only 
the self-bound gas and star particles within substructures are included in the visualization, therefore excluding the hot atmosphere of the cluster and both the 
diffuse intracluster stellar population and the stars associated with the BCG. The colours of the stars represent their age, using a red to light blue colour-table 
for decreasing age. Only a few of these galaxies, about 5-10 per cluster, still carries a self-bound hot gas halo, and those are located in the cluster peripheries. 
In the visualisation they appears as dark blue, extended halos, often with a comet-like shape, stretched away from the galaxies. 



portant already at distances of several virial radii (e.g. Dolag et al. 
2006). 

Using the "Zoomed Initial Conditions" (ZIC) technique (Tor- 
men et al. 1997), these regions were re-simulated with higher mass 
and force resolution by populating their Lagrangian volumes with 
a larger number of particles, while appropriately adding additional 
high-frequency modes drawn from the same power spectrum. To 
optimise the setup of the initial conditions, the high resolution re- 
gion was sampled with a 16'^ grid (for the more massive systems 
even with a 64^ grid), where only sub-cells are re-sampled at high 
resolution to allow for nearly arbitrary shapes of the high resolu- 
tion region. The exact shape of each final high-resolution region 



was iteratively obtained by repeatedly running dark-matter only 
simulations until the target objects were free of any contamina- 
tion by lower-resolution boundary particles out to 3-5 virial radii. 
The initial unperturbed particle distribution (before imprinting the 
Zeldovich displacements) was realized through a relaxed glass-like 
configuration (White 1996). 

Gas was then added to the high-resolution regions by split- 
ting each parent particle into a gas and a DM particle. The gas and 
the DM particles were displaced by half the original mean inter- 
particle distance, such that the centre-of-mass and the momentum 
of the original particle are conserved. As a further optimisation to 
save CPU time, the splitting of the DM particles in our largest sim- 
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ulations (i.e., those of the filamentary structure) was applied only 
to a subset of the high-resolution particles, selected by the follow- 
ing procedure. Instead of defining the high resolution region with 
the help of a grid, we used a surface around the centre of the high 
resolution region. For this purpose, we pixcliscd a fiducial sphere 
around the centre of the high resolution region using 4 x 16^ pix- 
els, making use of the HEALPIX tessellation of the fuU sky (G6rski 
et al. 1998), which has the advantage of offering an equal solid 
angle covered by each pixel. We then traced back in Lagrangian 
space all the particles in the DM-only run which at 2: = end up 
within 5 virial radii of the four massive clusters. All such particles 
were then projected onto the spherical coordinate system, where we 
recorded the maximum distance within each pixel of our HEALPIX 
representation. In a second step, all particles within the high reso- 
lution region were projected onto the spherical coordinate system 
and only those with distances smaller than the recorded maximvun 
distance of the corresponding pixel of our HEALPIX representation 
were then marked to define the region where to add gas. This op- 
timisation procedure reduced the number of gas particles by about 
40 per cent, without significantly reducing the size of the regions 
around the target objects where hydrodynamics is correctly com- 
puted. 

The final mass resolution of the gas particles in our sim- 
ulations is rrtgas = 2.4 x 10* M©. Thus, the massive clusters 
were resolved with between 2 x 10^' and 4 x 10^' particles. In 
all simulations, the gravitational softening length was kept fixed 
at £ = 42 kpc comoving Plummer-equivalent, and was switched to 
a physical softening length of e = 7 /t^^kpc at z = 5. 

The high resolution regions around the massive clusters are 
quite large and therefore also include other lower mass systems, 
which are still free of contamination from low-resolution particles. 
Hence, these systems can be included in our analysis. In this way, 
we end up with a total sample of 25 clusters having a mass of at 
least IO^^'Mq. The basic characteristics of these cluster, such as 
masses and radii at different overdensities, both for the DM-only 
runs and the non-radiative ovisc runs, are given in Table 2. 

To study resolution effects in the substructure mass distribu- 
tion in the DM runs, we also performed all these simulations at six 
times higher mass resolution, with the gravitational softenings de- 
creased accordingly by a factor 6^''^ ~ 1.8. In the following, we 
refer to these simulations as our high-resolution DM runs. 

Finally, In Appendix A we will show results of a resolution 
study of one single cluster, taken from a different set of cluster 
simulations (Borgani et al. 2006), which has been simulated at three 
different resolution. In the same Appendix we will also discuss the 
effect of adding gas particles in the initial conditions in such a way 
that there are 8 times fewer gas particles, but each having roughly 
the mass as that of the high-resolution DM particles. 



3 DETECTION OF SUBSTRUCTURES 

Substructures within halos are usually defined as locally over- 
dense, self-bound particle groups identified within a larger parent 
halo. In our analysis, the identification of these substructures is per- 
formed by applying the SUBFIND algorithm (Springel et al. 2001), 
which is defined in its original form only for dark matter only sim- 
ulations. In brief, as a first step, we employ a standard friends-of- 
friends (FoF) algorithm to identify the parent halos. In our anal- 
ysis we have used a FoF linking length of 0.16 times the mean 
DM particle separation. Note that this linking length is obtained 
when scaling the standard linking length of 0.2 by (Ac/fi)"^^^ 



according to the adopted cosmology and leads to groups having 
the overdensity characteristic of virialized objects predicted by the 
spherical collapse model (see Eke et al. 1996). The density of each 
particle within a FoF group is then estimated by adaptive kernel- 
interpolation, using the standard SPH approach with a certain num- 
ber of neighbouring particles. Using an excursion set approach 
where a global density threshold is progressively lowered, we find 
locally overdense regions within the resulting density field, which 
form a set of substructure candidates. The outer 'edge' of the sub- 
structure candidate is determined by a density contour that passes 
through a saddle point of the density field; here the substructure 
candidate joins onto the background structure. In a final step, all 
substructure candidates are subjected to a gravitational unbinding 
procedure where only the self-bound part is retained. If the number 
of bound particles left is larger than a prescribed minimum detec- 
tion threshold, we register the substructure as genuine subhalo. 

For full details on this substructure-finding algorithm as ap- 
plied to dark matter only simulations we refer the reader to the pa- 
per by Springel et al. (2001). In the following, we briefly describe 
the modifications to the algorithm that we implemented in order to 
make it applicable to simulations which also contain gas and star 
particles. 

Because the dark matter particles and the baryonic particles 
are in general distributed differently (especially in simulations with 
cooling), it is problematic to apply the FoF algorithm to all particles 
at once on an equal footing. This would not select groups that are 
bounded by a clearly defined density contour, and produces system- 
atic biases in the relative mass content of baryons and dark matter. 
We therefore apply the FoF only to the dark matter component of a 
simulation, using the same linking length that would be applied in a 
corresponding DM-only run. This selects halos that are bounded at 
approximately the same dark matter overdensity contour, indepen- 
dent of the presence or absence of baryonic physics. Each gas and 
star particle is then associated with its nearest DM particle, i.e. if 
this DM particle belongs to a FoF group, then the corresponding 
baryonic particle is also associated with that group. This effectively 
encloses all baryonic material as part of a FoF group that is con- 
tained in the original dark matter overdensity contour. 

We only keep FoF groups containing at least 32 DM particles 
for further analysis. Because a significant fraction of such small 
dark matter haloes are spurious, it is better to impose the detection 
limit just on the dark matter particle number, and not on the total 
number or total FoF group mass, because the latter would boost the 
number of spurious substructures detected in hydrodynamical runs. 
Also, this limit ensures that the halo number densities detected in 
DM-only and in hydrodynamic runs should only differ because of 
the effects of the baryonic component on the dark matter dynamics, 
and not because of a change in the group detection procedure. 

The second adjustment we made in SUBFIND is the procedure 
for density estimation in the presence of further components be- 
sides dark matter. We address this by first estimating the densities 
contributed by DM, gas particles and star particles at a given point 
separately for each of the components and then adding them up, 
i.e. the density carried by each of the three species is computed 
with the SPH kernel interpolation technique for neighbouring par- 
ticles of the same species only. Again, this is done to avoid possible 
systematic biases in the density estimates. If they are done with all 
particles at once, the different particle masses and spatial distribu- 
tions of the particle types can introduce comparatively large errors, 
for example when a few heavy dark matter particles dominate the 
sum over neighbours that mostly consist of light star particles. We 
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Cluster Sample 


N-4. 


a 


low mass 


141 ± 106 


-1.01 ±0.26 


low mass (high resolution) 


131 ± 32 


-0.97 ±0.15 


high mass 


136 ± 17 


-1.00 ±0.07 


high mass (high resolution) 


142 ± 22 


-0.99 ±0.03 



Table 3. The best-fit parameters for the cumulative subhalo mass function 
of Eq. (1). Reported here are the values for the low- and high-mass samples 
of DM-only simulations. Results are shown for the simulations performed 
both at the standard and at the high resolution. Quoted errors correspond 
to the rms scatter computed among the best-fitting values obtained for the 
individual clusters. 



note that we have carried out all the density computations with 32 
neighbours. 

Typically, the density contribution from the gas component is 
slightly smoother than that from the DM component (thanks to the 
gas pressure), whereas star particles generally trace a highly con- 
centrated density field, thus making it easier to detect substructures 
dominated by star particles. As a net result, we find that the result- 
ing total density field is a bit more noisy than the density field in 
pure DM runs. This led us to slightly modify SUBFIND's procedure 
to detect saddle points in the density field, and hence substructure 
candidates. In the default version of the algorithm, this is done by 
looking at the nearest two particles of a point that is added to the 
density field. We changed this to the nearest three points, finding 
that this leads to a significant improvement of the robustness of 
the identification of substructures in dissipative simulations while 
not changing the behaviour of the algorithm in pure DM or non- 
radiative runs. 

Finally, as a further refinement in SUBFESTD, we take into ac- 
count the internal thermal energy of the gas particles in the grav- 
itational unbinding procedure. After unbinding, we keep substruc- 
tures in our final list of subhalos if they contain a minimum of 20 
dark matter and star particles. We ignore gas particles in this detec- 
tion threshold in order to avoid detection of spurious gravitationally 
bound gas clumps. In principle, a cleaner mass threshold for detec- 
tion would be obtained by also ignoring the star particle count for 
this validation, but including them allows us to identify as genuine 
self-bound substructures also very low mass systems that are dom- 
inated by the stellar component. Since stars form at the bottom of 
the potential wells and are highly concentrated relative to the other 
two components, we find this does not produce a significant com- 
ponent of spurious subhalos. 

For each of the FoF groups that correspond to a cluster, we 
obtain in this way a catalogue of gravitationally botmd subhalos 
that in general contain a mixture of dark matter, star and gas par- 
ticles. The centre of each subhalo is identified with the minimum 
of the gravitational potential occuring among the member parti- 
cles. Similarly, for the cluster as a whole, we determine a suitable 
centre as the position of the particle that has the minimum grav- 
itational potential. Around this point, we calculate the virial ra- 
dius and mass with the spherical-overdensity approach, using the 
over-density predicted by the generalized spherical top-hat collapse 
model (e.g. Eke et al. 1996). We define as subhalos of a cluster all 
the substructures which are identified within its spherical virial ra- 
dius. Note that this can sometimes include subhalos that belong to 
a FoF group different from that of the main halo, due to the aspher- 
ical shapes of the FoF groups. Likewise, not all of the subhalos in 
the cluster's FoF group are necessarily inside its virial radius. 



4 SUBHALOS IN GALAXY CLUSTERS 
4.1 DM-only simulations 

As a first step, we characterize the mass distribution of subhalos 
for the DM runs and check the effect of the mass of the parent 
halo and of resolution. For this purpose, we split our sample of 
simulated clusters into a high mass sample, containing 8 halos with 
virial masses above 10^' M©, and a low mass sample, containing 
13 halos with masses ranging between ~ 1.0 x IO^^Mq and « 
3.0 X 1O"M0 (see Table 1). 

Figure 3 shows the cumulative subhalo mass function, normal- 
ized to the virial mass of the clusters (i.e. subhalo masses are mea- 
sured in units of the virial mass of the parent halo). In the left panel 
and in the central panels, the Unes are for individual clusters, with 
grey and black lines showing the low-mass and high-mass samples, 
respectively. The scatter among different clusters is substantial for 
substructure counts below ^ 100, which only includes the most 
massive structures. However, once a few hundred substructures are 
considered, the abundance per unit virial mass becomes quite uni- 
form. There appears to be no significant difference in the scatter of 
the high-mass end of the subhalo mass functions between our low- 
and high-mass cluster samples. 

In the middle panel, we show the same results as in the left 
panel, but now for 6 times higher mass resolution. This allows us 
to coimt down to subhalos that are 6 times smaller in mass, yield- 
ing significantly higher total subhalo counts. There is no clear dif- 
ference in the cluster-to-cluster scatter relative to the lower resolu- 
tion simulations, suggesting that the scatter is not due to numerical 
noise but rather reflects the genuine variations in the abundance of 
massive subhalos in different halos. 

This is corroborated by the right panel of Fig. 3, where we 
compare the average mass functions for the two subsets of low- 
mass and high-mass clusters, and at the two numerical resolutions. 
The thick grey line marks a power law with slope of —1, for com- 
parison. Clearly, the average slope of the subhalo mass function of 
both cluster samples and at both resolutions agree perfectly. 

This result is good in agreement with previous findings (e.g. 
Gao et al. 2004; De Lucia et al. 2004). However, unlike Gao et al. 
(2004), we here do not find a significant trend of the amplitude 
of the subhalo mass function as a function of halo mass, a results 
which is in line with the findings by De Lucia et al. (2004). How- 
ever, it may simply be the relatively small range of halo masses 
probed by our simulations that prevents us from seeing the weak 
mass trend, despite the good resolution reached in our most mas- 
sive systems. We note that in the high-resolution versions of the 
high-mass clusters, each main halo is resolved with more than 10 
million dark matter particles within i?vir, allowing several thousand 
subhalos to be identified within each main halo. 

A good quantitative fit to our measured subhalo mass func- 
tions can be obtained with a power-law of the form 



(1) 



We obtain values of the slope a and of the normalization N-4 for 
each cluster by performing a fit to the numerical data above a lower 
limit of 100 DM particles (to be insensitive to the resolution limit) 
and for msub/A/vir < 0.01 (in order to reduce sensitivity to the 
scatter in the high-mass end of the subhalo mass function). The re- 
sulting fitting parameters arc reported in Table 3, along with mean 
values over our cluster sample and the errors due to the object- 
by-object scatter. The results of this table confirm that at high res- 
olution we obtain rather stable results for both the slope and the 
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Figure 3. The cumulative subhalo mass functions for the DM runs of our clusters, normalized to the virial mass of each cluster. Left panel: comparison 
between the 13 low mass halos (in grey) and the 8 high mass halos (in black). Central panel: the same as in the left panel, but for the runs with six times higher 
mass resolution. Right panel: the average cumulative subhalo mass function for the low (in grey) and high (in black) mass sample from the DM-only runs. The 
dashed Hnes are for the simulations at the standard resolution while the solid lines are for the higher resolution runs. The coiTesponding best-fit values of the 
power-law fit to each of these curves are given in Table 3. For reference, the grey thick line marks a power-law with slope —1. 
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Figure 4. The cumulative total mass functions (left and central panel) and stellar mass functions (right panel) for the subhalos identified in the g51 cluster at 
z = 0, as obtained for different runs. Left panel: comparison of the total subhalo mass functions of the DM run with the non-radiative runs corresponding 
to different schemes for the artificial viscosity (solid lines). The dashed lines are obtained by correcting the mass of each sub-halo in the non-radiative runs 
by increasing it by the cosmological baryon fraction. Central panel: the same as the left panel, but for the different radiative runs. Right panel: comparison 
between the stellar mass functions obtained for the different radiative runs. As in Figure 3, the thick grey line marks a power-law with slope — 1. 



normalisation of the subhalo mass functions. Furtheimore, we find 
no significant differences between the mass functions of high and 
low mass clusters, both in slope and in normalization. 

4.2 Hydrodynamical runs 

We now turn to an analysis of the subhalo mass functions in our 
hydrodynamical simulations. As described in Section 3, we have 
modified the subhalo detection algorithm SUBFIND such that it can 
properly take into account the presence of gas and star particles, be- 
sides the dark matter particles. Our modifications have been guided 
by the desire to avoid possible systematic biases due to the pres- 
ence of multiple particle components, such that the final subhalo 
catalogues can be directly compared with those of a DM-only sim- 
ulation, except that the subhalos can now also contain gas and star 
particles. 

Figure 4 shows a comparison of the cumulative subhalo mass 



functions for the different runs we carried out for the high-mass 
gSl.a cluster. The left panel compares the total subhalo mass func- 
tion for the DM-only run with three non-radiative hydrodynamic 
runs (solid lines) that used different treatments for the numerical 
viscosity (see Section 2.1 for details). Clearly, introducing non- 
radiative hydrodynamics causes a decrease of the total subhalo 
mass function by a sizable amount, and this is almost independent 
of the detailed scheme used for the artificial viscosity. If interpreted 
in terms of a shift in mass, the difference between the subhalo mass 
in the DM-only and in the gas runs can be explained mostly by the 
high efficiency of gas removal during the infall of subhalos into the 
cluster. To demonstrate this, we increased the mass of each sub- 
halo identified in the non-radiative run by an amount correspond- 
ing to the cosmic baryon fraction. If gas is completely stripped and 
this represents the only reason for the lower mass function, then 
this correction should provide a mass function consistent with that 
of the corresponding DM run. However, the stripped mass can not 
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Figure 5. Left panel: the concentration parameter, obtained by fitting a NFW profile (Navarro et al. 1996) to the dark matter density profiles in the hydrody- 
namical simulations, as a function of J\/200- Here we used all halos with masses above 10^^ Mq identified in all the simulation snapshots at redshift above 
z = 0.5, thereby scaling the obtained concentration parameter by (1 + z) to compensate for the evolution. The different lines mark the results for the dark 
matter run (dm), the non-radiative run (ovisc) and the run with cooling and star formation (csf). The dashed line gives the fit proposed in Dolag et al. (2004). 
The right panel shows the half mass radius as a function of the sub-halo mass for all substructures in all the clusters at redshift 2 = 0. The different lines 
distinguish the different simulations as before. In both panels the eiTorbars (shown only for the dm run) correspond to the rms scatter among all the clusters 
within each mass bin. 




Figure 6. Evolution of the subhalo mass functions for the different runs, averaged over the 8 clusters of the high mass sample. The two left and the two central 
panels show the mass functions at z = 0, 0.3, 0.5 and 1 with sohd, dashed-dotted, dashed and dotted lines, respectively. Upper left, upper central and lower 
left panels are for total mass function for the DM, non-radiative ovisc and radiative ci/runs, respectively. The lower central panel is for the stellar mass function 
of the radiative csf runs. The two right panels compare the total mass function for DM, ovisc and csf runs at z = (upper right panel) and at z = 1 (lower 
light panel). The thick grey hues in all panels mark the slope —1. 
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completely compensate for the mass difference between the dark 
matter and the non-radiative runs. Besides removing gas, also their 
orbits are changed due to the effect of gas pressure (see also Tor- 
men et al. 2004; Puchwein et al. 2005; Saro et al. 2008), also with 
some indications that the remaining dark matter sub-halos might 
become easier to disrupt. 

The central panel of Figure 4 shows the cvmiulative subhalo 
mass functions for the same cluster, but here comparing the DM 
run with the different radiative runs, which differ in the efficiency 
of the kinetic feedback (i.e. csfnw, csf, and csfsw), or in the pres- 
ence of thermal conduction (i.e., csfc). The bulk of gas cooling and 
subsequent star formation within subhalos takes place before the 
galaxies are falling into the high-pressure environment of the ICM. 
Cooling has the effect of concentrating cold gas and therefore also 
collisionless stars forming out of the gas at the centres of these sub- 
halos. This increases the concentration of the mass distribution, and 
protects part of the baryonic mass from ram-pressure stripping. As 
a result, the subhalo mass functions in these radiative runs increase 
when compared to the non-radiative cases and become quite sim- 
ilar to those of the DM-only runs, with only a marginal sensitivity 
to the adopted intensity of the feedback. Only in the run without 
any kinetic feedback {csfnw), the resulting total subhalo masses are 
even slightly larger than for the DM runs. This is probably the con- 
sequence of the adiabatic contraction of the DM halos, in response 
to the rather strong cooUng taking place when galactic winds are 
turned off (e.g., Gnedin et al. 2004). 

It is worth mentioning that the change of the disruption of 
sub-halos in the hydrodynamical simulations is due to a complex 
interplay between different processes. Firstly, we find that the pres- 
ence of baryons makes the underlying dark matter distribution more 
concentrated (see also Rasia et al. 2004; Lin et al. 2006). A further 
increase in the concentration of the dark matter profiles is also con- 
tributed by including radiative cooling. These findings are reported 
in the left panel of Figure 5, which shows the concentration param- 
eter of dark matter profiles as a function of halo mass for all main 
halos. The results for the dark matter {dm) runs are in good agree- 
ment with the fit originally presented by Dolag et al. (2004). The 
concentrations in the non-radiative runs and in the runs with cool- 
ing and star-formation increase on average by 3-8% and by 10-25% 
with respect to the dm runs, respectively. The expectation based 
on this result is then that the halos in the non-radiative runs are 
more difficult to destroy. However, as soon as the gas component is 
stripped from the sub-halos, the latter react to this mass loss with 
a slight expansion. In the non-radiative runs {ovisc) this expansion 
even slightly over-compensates the increase of the concentration. 
This is demonstrated by the right panel of Figure 5, which shows 
the half-mass radii of all substructures found in the main haloes at 
a = as a function of their total mass. We find that this radius in 
the non-radiative runs {ovisc) is ~ 2% larger than in the pure dm 
runs. However, in the runs with cooling and star formation {csf), 
the sub-halos remain more concentrated and we find the half-mass 
radius to be « 5% smaller than in the dm run. This increase in 
halo concentration is sufficient to compensate for the change of or- 
bits due to the drag induced by the ram pressure (see Tormen et al. 
2004; Puchwein et al. 2005; Saro et al. 2008). As a consequence, 
sub-haloes survive for a longer time in the csfmws, whereas these 
effects are nearly compensating in the non-radiative runs with a 
possible tendency for sub-haloes to be more easily disrupted when 
compared to the dm only runs. 

For the radiative runs, we note that a small but sizable popu- 
lation of small subhalos exist that are dominated by star particles 
and whose DM component was lost during the subhalo's evolution. 



These surviving compact cores of star particles are still recognised 
by SUBITND as self-bound structures, and give rise to the low- 
mass extension of the subhalo mass functions. The flattening of the 
mass function in these regimes shows that the number of DM-poor 
galaxies is affected by the finite numerical resolution. 

The right panel of Figure 4 shows the subhalo mass functions 
of just the stellar component, for the different radiative runs. As ex- 
pected, the effect of changing the feedback strength is now clearly 
visible and can lead to nearly a factor of ten in stellar mass between 
simulations without kinetic feedback {csfnw) and with strong ki- 
netic feedback {csfsw, see also Saro et al. 2006). This comparison 
highlights the existing interplay between the feedback, which reg- 
ulates the star-formation within subhalos and its progenitor halos, 
and the stripping due to the highly pressurised environment of the 
intra-cluster medium. Although thermal conduction is expected to 
alter the efficiency of gas stripping (e.g., Price 2007), its low effi- 
ciency in the low-temperature subhalos results only in a marginal 
impact on the resulting stellar mass function. In Appendix A, we 
will further examine the numerical convergence of the stellar mass 
function of the cluster galaxies. 

In Figure 6, we consider the time evolution of the subhalo 
mass functions by showing averages for the main progenitors of 
the eight high-mass clusters at different redshifts. For the DM-only 
simulations we do not see any significant evolution over the con- 
sidered redshift range. This is not directly in contrast to findings 
by Gao et al. (2004), as they define the vuial mass with respect to 
200 times the critical density. If we stick to the same definition of 
virial mass and radius we also get a evolutionary trend in the same 
direction than found in Gao et al. (2004). The analysis of the aver- 
age subhalo mass function over the full sample of massive clusters 
confirms the effect of the gas stripping in the non-radiative runs, al- 
ready discussed for the gSl cluster at 2 = 0, and the counteracting 
effects by cooling and star-formation in the radiative runs. 

Furthermore, the rightmost upper and lower panels show the 
results at 2 = and 2=1, respectively. At intermediate sub- 
halo mass (between 10"Mc5 and IO^^Mq), the effects of coohng 
and star formation are that of slightly increasing the subhalo mass 
function compared with the pure DM run. This is especially true at 
high redshift, when stripping is less efficient. There are indications 
that, at the high mass end, stripping can be still efficient enough 
to suppress the subhalo mass function, also in the runs with cool- 
ing and star formation. Again, a measurable effect is seen from 
star dominated substructures for the low mass subhalos, indicating 
that clumps of star particles, being more compact than DM sub- 
halos, are more resistant against tidal destruction. However, even 
the radiative simulations do not predict that a significant fraction of 
galaxies without dark matter halos survive down to 2; = 0. 

Figure 7 shows the evolution of the parameters of the power- 
law fit to the total subhalo mass function, where we have as usual 
restricted the fitting range to halos with at least 100 dark matter 
particles and msub/A/vir < 0.01, which gives a reasonable fit over 
the whole considered redshift range. To guarantee robust fits, we 
here excluded those halos for which there are less than 20 data- 
points in the estimate of the cumulative mass function over this 
range. The slope of the total subhalo mass function does not show a 
significant evolution in all the runs. As for the normalisation (right 
panel of Fig. 7), its value for the non-radiative {ovisc) gas run is 
always below that of the DM and of the radiative csfgss runs. There 
is also an indication that the radiative run at high redshift has a 
normalization slightly higher than the DM run. The stripping of 
the gas in the non-radiative run is relatively independent of the halo 
mass, thus leading to a slope similar to that of the DM case. We note 
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Figure 7. The evolution of the parameters defining the power-law shape of the subhalo mass function (see equation 1) for the progenitors of the 8 clusters 
belonging to the high-mass subsample. All the 52 outputs of the simulations out to z = 1.4 have been used to produce the continuous lines. Different 
colours/symbols are for the DM runs, for the {ovisc) non-radiative run and for the {csf) radiative runs. The left and the right panels show the evolution of slope 
and normalisation, respectively. The error-bars mark the rms scatter computed over the 8 clusters of the high-mass subsample. For reasons of clarity, they are 
shown along with the symbols only at z = 0, 0.3, 0.5 and 1.0, and are slightly displaced in the horizontal direction. 



that there is a tentative indication that the condensation of baryonic 
matter at the centre of the substructures in radiative runs leads to 
a small systematic trend with redshift, leading to a slightly flatter 
subhalo mass function at early redshift in the c.s/runs. 



4.3 The composition of subhalos 

In DM-only simulations, it is known that tidal stripping induces 
strong radial dependences in the properties of substructures (e.g. 
Gao et al. 2004; De Lucia et al. 2004). For example, the substruc- 
ture abundance is antibiased relative to the mass distribution of 
the cluster, and the average concentration of subhalos increases 
towards the cluster centre. The efficiency of gas stripping from 
substructures also depends strongly on cluster-centric distance (see 
Gunn & Gott 1972; McCarthy et al. 2008), hence we expect radial 
variations in the relative content of baryons in subhalos, which we 
analyse in this section. 

The problem is made especially interesting by the complex 
non-linear interplay between the condensation of baryons due to 
radiative cooling, the gas-dynamical stripping processes, and the 
resulting orbital evolution of affected subhalos. Indeed, while gas 
stripping in the high-pressure ICM makes subhalos more fragile, 
the formation of a compact core of cooled gas and stars makes 
subhalos more resistant to disruption. As a result, we expect that 
the relative stellar content of a subhalo and its capability to retain a 
gaseous halo will depend both on the details of the physics included 
in the simulation and on its cluster-centric distance. 

In order to investigate this in more detail, we show in Figure 8 
the radial dependence of the fraction of subhalo masses contributed 
by stars for the radiative runs including the reference kinetic feed- 
back {csf), at three different epochs {z = 0, 1, 2 from left to right). 
In each panel, the different lines indicate different limits for the 
total mass of the selected subhalos. Quite clearly, the stellar com- 
ponent becomes progressively more important in subhalos found 
at a smaller cluster-centric distance. This is consistent with the ex- 
pectation that the strong tidal interactions in the cluster central re- 
gions are efficient in removing the dark matter component of the 
subhalos, while the gravitationally more tightly bound clumps of 
star particles are more resistant. Fig. 8 also shows that the radial 



dependence of the subhalo mass fraction in stars extends well be- 
yond the clusters' virial radius i?vir, thus indicating that the cluster 
environment affects the properties of the galaxy stellar population 
over a region much larger than the virial one (see also Saro et al. 
2006). Note that a substantial fraction of these subhalos detected at 
r > Rvir may have already crossed the cluster virial region in ear- 
lier passages and constitute the so-called "backsplash population" 
(see Gill et al. 2005; Ludlow et al. 2008). Their orbits have brought 
them into the outer parts again, with a fraction of them eventually 
destined to escape the cluster. 

Towards the centre of the clusters, the star fraction within 
subhalos reaches values higher than the cosmic baryon fraction. 
At first glance, this result lends support to the assumption made 
by many semi-analytical models (SAMs) of galaxy formation that 
galaxies can preserve their identity for a while even after their host 
DM-halos is destroyed by the tidal interaction within clusters (Gao 
et al. 2004). However, in semi-analytical models (like De Lucia & 
Blaizot 2007) based on DM simulations of comparable resolution 
(like the Millennium Run, Springel et al. 2005) about half of the 
galaxies within the virial radius of the parent cluster are expected 
to have lost their DM halos at z — 0, while in our hydrodynamical 
simulations we find that only ~ 20 per cent of the galaxies have 
a star component which exceeds the mean cosmic baryon fraction 
(i.e. m*^^, > 0.13 msub). while only 2 per cent of the galaxies have 
their mass dominated by star particles (i.e., m*^^, > 0.5maub). This 
implies that only a small fraction of our star dominated systems can 
be classified as orphan galaxies, which are defined as those galaxies 
not having a counterpart in the corresponding dark matter only run. 
A precise comparison requires however a carefully matched mag- 
nitude selection, which is outside the scope of the present paper. 
Future, yet higher resolution hydrodynamical simulations should 
in principle allow an accurate calibration of the SAM assumptions 
about the relative survival time of satellite galaxies once their par- 
ent dark matter halo cannot be tracked any more. 

A related question concerns the ability of infalling subhalos to 
retain a gravitationally bound component of diffuse gas during the 
cluster evolution. We show in Figure 9 the cluster-centric depen- 
dence of the subhalo mass fraction associated with gravitationally 
bound gas (in units of the cosmic baryon fraction) from the set of 
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eight massive clusters, at 2; = 0, 1 and 2. In each panel, we compare 
results for the ovisc non-radiative run with those for the radiative 
cj/run. Furthermore, the solid lines indicate the radial dependence 
of the average subhalo gas mass fraction, which is computed by in- 
cluding also those subhalos which do not have any gravitationally 
bound gas component. We indicate with the thin lines the Poisson 
uncertainties in the estimate of the average gas mass fraction, by 
assuming them to be proportional to the square-root of the number 
of galaxies in each radial bin. At 2 = we note that, whenever 
a subhalo retains a gas component, the associated mass fraction is 
higher in the non-radiative run than in the radiative one. This is 
related to the fact that a significant fraction of this gas is converted 
into stars in the radiative case. 

On the other hand, the central concentration of star particles 
in subhalos in the radiative runs makes them more resistant against 
ram-pressure stripping. As a consequence, a larger number of sub- 
halos preserve a gas component, thus explaining the higher average 
gas fraction at all radii. Within the virial radius, we find at 2 = 
that about ~ 1 per cent of the resolved subhalos still have some 
self-bound gas, both in the ovisc and in the ci/runs, with a marginal 
indication for a smaller amount of gas in the «/ clusters, due to its 
conversion into stars. As expected, at high redshift ram-pressure 
stripping has not yet been as effective, thus explaining the larger 
fraction of self-bound gas, which is ~ 4 per cent and ~ 6 per 
cent for the ovisc and the c.s/runs, respectively. Again, we note that 
the radial dependence of the gas fraction extends well beyond the 
virial radius, thus confirming that the overdense cluster environ- 
ment affects the structure and evolution of subhalos over a fairly 
large region surrounding the clusters (see also Dolag et al. 2006). 
Also note that at very high redshift the collapse of the gas due to 
efficient cooling inside the subhalos can even lead to gas fractions 
larger than the cosmic values in the radiative runs. 

An important question for hydrodynamical simulations of cos- 
mic structure formation is how faithfully they are able to describe 
the gas-dynamical processes which determine the stripping of the 
gas in substructures moving in a high density environment. In these 
processes, hydrodynamical fluid instabilities that are difficult to re- 
solve can play an important role. In general, the faithfulness of the 
simulation depends both on the hydrodynamical scheme adopted 
(i.e., Eulerian or Lagrangian, implementation of gas viscosity, etc.) 
and on the numerical resolution achieved (e.g., Dolag et al. 2005; 
Sijacki & Springel 2006; Agertz et al. 2007; lapichino & Niemeyer 
2008). Although a detailed study of these numerical aspects is be- 
yond the scope of this paper, it is interesting to examine the impact 
of the adopted artificial viscosity parameterization on the gas strip- 
ping from subhalos. 

Figure 10 shows the subhalo gas mass fraction as a function of 
cluster-centric distance for the four clusters of our high-mass sam- 
ple, for which runs with different artificial viscosities have been 
carried out (i.e., simulations gl.a, gS.a, gSl.a and g72.a\ see Ta- 
ble 2. Similarly to the results obtained by Dolag et al. (2005) for 
the predicted ICM turbulence, the ovisc and svisc implementations 
give quite similar results, whereas the Ivisc low-viscosity scheme 
leads to smaller gas-mass fractions inside the halos and slightly 
higher gas-mass fractions outside the virial radius. This is a clear 
indication that numerical viscosity has a measurable effect on the 
stripping of gas from the subhalos, as the latter depends on the en- 
vironment. 

The different behaviour of the lower-viscosity scheme can be 
understood by recalling that this approach provides less efficient 
viscous stripping, which explains the larger gas fraction associated 
with subhalos found in the cluster outskirts. On the other hand. 
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Figure 8. The evolution of the radial profiles of the stellar mass fraction 
within subhalos for the radiative csf run. Each curve is obtained by aver- 
aging over the 8 clusters belonging to the high-mass sample. Only subha- 
los with masses ^ 2 X 10^^ Mq have been included in this analysis. The 
horizontal dashed line marks the cosmic baryon fraction assumed in our 
simulations. 

as sub-halos cross the virial radius and enter a progressively more 
pressurised environment, they should experience hydrodynamical 
instabilities in the shear-flows around the subhalo's gas, such as the 
Kelvin-Helmholtz instability, which may be poorly represented in 
SPH (Agertz et al. 2007). While the low-viscosity scheme does not 
improve the SPH description to a point where these instabilities are 
captured satisfactorily with SPH, it still improves the description of 
shear flows, so that the subhalo gas can be dissolved and dispersed 
in the hotter atmosphere over a relatively short time-scale. 



5 CONCLUSIONS 

In this study, we presented an analysis of the basic substruc- 
ture properties in high-resolution hydrodynamical simulations of 
galaxy clusters carried out with the TreePM-SPH code GADGET. 
We used a fairly large set of 25 clusters, with virial masses in the 
range 1 — 33 x 10^'' /i^^Mq, out of which eight clusters have 
masses above 10^^ K~^Mq. The identification of the substructures 
has been carried out with the SUBFIND algorithm (Springel et al. 
2001), which we suitably modified to account for the presence of 
gas and star particles in hydrodynamical simulations. The primary 
aim of our analysis was to quantify the impact that gas physics has 
on the properties of subhalos, especially on their abundance. Note 
that previous work has so far almost exclusively analyzed substruc- 
tures in dark matter-only simulations, so that our study is one of the 
first attempts to directly compare the DM-only results to those of 
hydrodynamical simulations that also account for baryonic physics. 

We analyzed non-radiative hydrodynamical simulations as 
well as radiative simulations, including star formation and different 
efficiencies for energy feedback associated with galactic outflows. 
Also, we examined the sensitivity of our results with respect to dif- 
ferent parameterizations for the artificial viscosity needed in SPH. 
The main results of our analysis can be summarized as follows. 

• Consistent with previous studies (Springel et al. 2001; De Lu- 
cia et al. 2004), we find that the subhalo mass function in DM- 
only runs converges well for different numerical resolution over 
the mass range where substructures are resolved with at least ~ 30 
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Figure 9. The radial dependence of the fraction of subhalo mass in gas. in units of the cosmic baryon fraction, for the eight clusters belonging to our high-mass 
sample. Black and blue symbols are for the subhalos retaining a gravitationally bound gas component for the non-radiative ovisc and for the radiative c.s/runs, 
respectively. The heavy continuous lines show the radial dependence of the gas fraction after averaging over all subhalos. including those not retaining any gas. 
while the light lines indicate the uncertainties arising from the small number statistics of the involved subhalos. Left, central and right panels are for z = 0, 1 
and 2, respectively. 




Figure 10. The same as in Figure 9, but comparing results for the different implementation of ailificial viscosity. Results are shown for the four clusters of the 
high-mass sample for which runs with the different viscosity schemes have been carried out. 



gravitationally bound DM particles. We extend this result to show 
that the same convergence also holds well for the total subhalo mass 
function in hydrodynamical runs, irrespective of whether cooling 
and star formation are included in the simulations. 

• In non-radiative hydrodynamical simulations, the presence of 
a gas component leads to a reduction of the total subhalo mass 
function with respect to DM-only runs. On average, the reduction 
in mass of the substructures is slightly larger than expected from 
the complete stripping of the baryonic component. This is due to a 
combination of the modification of the subhalo orbits arising from 
the pressure force exerted by the intra-cluster medium. Further- 
more, we see indications for an increased susceptibility of the sur- 
viving DM halos to tidal disruption. 

• In general we find that only a very small fraction of subhalos, 
of order of one percent, within the cluster virial radii maintain a 
self-bound atmosphere of hot gas. This indicates a high efficiency 
for the stripping of the hot gas component in our simulations. The 
radial dependence of the mean gas fraction within subhalos indi- 
cates that gas stripping starts already at large cluster-centric dis- 
tances, beyond the virial radius. Using a scheme for reduced arti- 



ficial viscosity has the effect of reducing viscous stripping in the 
cluster periphery, while increasing the stripping within the virial- 
ized high-density region. 

• In radiative runs that include star formation, some of the 
baryons are protected from stripping due to their concentration at 
the centres of the subhalos. In this case, substructures have masses 
which are comparable to, or are slightly larger than in the DM-only 
runs. The stellar mass fraction is also found to strongly increase for 
subhalos close to the centre. This confirms that compact clumps 
of star particles are indeed more resistant against tidal destruction 
than DM subhalos. Despite this, we find that only a small fraction of 
subhalos are dominated by their stellar component. As a result, the 
number of star-dominated galaxies is smaller than expected when, 
like usually assumed in semi-analytical models of galaxy forma- 
tion, they are allowed to survive for a dynamical friction time after 
the disruption of their DM halo. This suggests that our hydrody- 
namical simulations at their current resolution may not be able yet 
to accurately follow the survival of galaxies within clusters after 
the destruction of their DM subhalos. 

Our study shows that the technical improvements we imple- 
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mented in SUBFIND allow a reliable identification of substructures 
in hydrodynamical simulations that include radiative cooling and 
star formation. This allows a direct identification of satellite galax- 
ies not only in clusters of galaxies, but also in future high-resolution 
hydrodynamical simulations of galaxy-sized halos. This is essen- 
tial to make direct contact with observational data, and also allows 
studies of the structural impact of baryonic physics on substruc- 
ture. As our result show here, gas physics does alter substructure 
statistics in interesting ways, but the effects depend on physical 
processes such as radiative cooling and feedback. This also means 
that the approximate self-similarity of substructure properties that 
has been found in DM-only simulations of clusters and galaxy-size 
objects (Moore et al. 1999) is not expected to hold nearly as well 
in simulations with radiative cooling. For example, in galaxy-size 
halos, many subhalos will be so small that they do not experience 
atomic cooling; as a result their behaviour should be close to our 
non-radiative results, and we would here expect a reduction of their 
abundance relative to the DM-only case. 

In future work, it will be important to further improve the nu- 
merical description and resolution of our simulations, in order to be 
able to more accurately follow in particular the star-dominated sub- 
structures (i.e., galaxies) within galaxy clusters. This will then also 
provide important input for the semi-analytical models of galaxy 
formation, allowing a check and calibration of their dynamical fric- 
tion and gas stripping prescriptions. 
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APPENDIX A: NUMERICAL STABILITY 

This Appendix is devoted to a discussion of the numerical stabil- 
ity of the results presented for the properties of the substructures 
in simulated clusters. In the first part, we will discuss the robust- 
ness of the results against the choice of the subhalo finding algo- 
rithm. In the second part we will discuss the stability of the sub- 
halos mass functions against numerical resolution, choice of the 
minimum SPH smoothing length and scheme to add gas particles 
to initial conditions. 

Al Substructure detection 

In Section 3 we described our modification of the SUBFIND algo- 
rithm (Springel et al. 2001) to allow a clear detection of substruc- 
tures with various compositions in hydrodynamical simulations in- 
cluding cooling and star formation. To verify our method, we rean- 
alyzed some of the clusters studied in Murante et al. (2007) with the 
new method and compared the identified stellar components of the 
individual galaxies. Note that the identification of galaxies in Mu- 
rante et al. (2007) was based on a completely independent method, 
namely SKID (Stadel 2001). Note that SKID works best if applied 
only on the distribution of star particles. Therefore, for the sake of 
comparison, we excluded from our SUBFIND galaxies those having 
no stellar component. 

Figure Al shows a comparison for the R3 run of the positions, 
individual masses and stellar subhalo mass functions. We find an 
excellent agreement between the stellar masses inferred from both 
algorithms, in fact, the majority are identical for both algorithms. 
Only for the very loosely bound star particles at the periphery of 
the substructures some differences occur. The middle panel of the 
figure shows a comparison between the stellar masses of the galax- 
ies identified by the two algorithms. Note that the horizontal and 
the vertical branches at low masses indicate the galaxies which are 
identified by one of the two algorithms. However, it is worth men- 
tioning that there is no indication of any systematic offset between 
the results, and most of the scatter arises because sometimes merg- 
ing substructures get split into two parts only by one of the two al- 
gorithms. In this respect, we find that SUBFIND is marginally more 
efficient in separating two merging substructures. As a general re- 
sult of our comparison, we conclude that our modifications of the 
SUBFIND algorithm are able to reliably identify the stellar parts of 
subhalos even though we apply the substructure finder to all the 
particles in the simulations. 

A2 Resolution and convergence 

To test our results we performed several DM-only runs by a chang- 
ing a number of parameters defining the numerical setup. Among 
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Figure Al. Comparison of the galaxies detected by SKID and by SUBFIND for the R3 run at z = 0. The left panel shows the position of the individual 
substructures, with different symbols/colours for galaxies identified with the two algorithms, as specified in the labels. Red diamonds and black crosses indicate 
galaxies identified by SUBFIND and by SKID, respectively, when the positions provided by the two algorithms differs by less than the half mass radius of 
the galaxy in the SUBFIND identification. Blue diamonds and green crosses are for the galaxies identified by SUBFIND and by SKID, respectively, when 
they do not have a counterpart identified by the other algorithm. In the results from SUBFIND, we do not include identified substructures which do not have a 
stellar component. The middle panel shows a one-to-one comparison between the stellar masses assigned by the two algorithms to the identified galaxies. The 
symbols lining up veilically/horizontally are those identified by one algorithm, but with no match in the galaxy sample from the other algorithm. The dashed 
lines mark the stellar mass limit applied to the SKID detection. The right panel shows the cumulative stellar subhalo mass function from the two identification 
schemes. We also plot the mass functions for the galaxies identified from one algorithm and having a match in the galaxy sample derived from the other 
algorithm. 



the tests carried out, we mention the following: (a) use only the 
Tree-based gravity solver and compare with the Tree-PM scheme; 
(b) increase the accuracy of time integration and force computa- 
tion far beyond the usual values assumed in our runs (c) change 
from z = 5 to z = 2 the redshift at which the gravitational soften- 
ing switches from physical to comoving units. In all these tests we 
verified that the subhalo mass functions do not show appreciable 
variations down to the mass corresponding to 20 self-bound DM 
particles. 

As already discussed in Section 4.1, increasing resolution in 
the DM runs produces numerically stable results in the subhalo 
mass function, at least as long as only structures with at least 20 
gravitationally bound DM particles are selected. The situation is 
expected to be different when considering instead hydrodynamical 
simulations, which include physical processes that are expected to 
be much more resolution dependent. The clearest example is ra- 
diative cooling, whose efficiency is known to depend on the mass 
resolution in cosmological simulations. In a similar fashion, hydro- 
dynamical processes, such as ram-pressure and viscous stripping, 
depend quite sensitively on the resolution adopted. 

In order to assess the effect of resolution on the subhalo mass 
function, we carried out several simulations of a moderately poor 
clusters, having a virial mass Mvir ^ 2.3 X W^Mq. The set-up of 
the initial conditions of this cluster have been described by Borgani 
et al. (2006) (it corresponds to the CL4 cluster of that paper). The 
cosmological model assumed in this case is the same as that for the 
sample of clusters analysed in this paper, except that a lower nor- 
malisation (78 = 0.8 of the power spectrum was assumed (Borgani 
et al. 2004). Radiative simulations of this cluster have been carried 
out by varying the effective resolution in three different ways: 

• Increase the mass resolution by a factor 4 and by a factor 
15 with respect to a reference resolution. The corresponding force 
resolutions are decreased by a factor 4^^^ and IS'^'^^, respectively 
(runs Rl, R2 and R3 in Table Al). 



• Change the criterion to assign gas particles in the initial con- 
ditions. Instead of using the same number of gas and DM particles, 
with a mass ratio reflecting the cosmic baryon fraction, we assigned 
one gas particle to every group of eight DM particles. In this way, 
the mass of each parent DM particle is decreased by an amount 
such to reproduce the cosmic baryon fraction. Given the standard 
values of the cosmic baryon fraction, this procedure gives compa- 
rable masses to DM and gas particles, thus reducing a possible spu- 
rious numerical heating of the gas particles induced by two-body 
encounters with DM particles. This scheme of adding gas to ini- 
tial conditions can be seen as a way of increasing the accuracy of 
the gravity part (i.e. increasing the number of DM particles) for a 
fixed resolution in the hydrodynamics. This prescription was origi- 
nally suggested by Steinmetz & White (1997) as a computationally 
efficient way of suppressing two-body heating in cooling flows. 

• Use different values for the minimum SPH kernel smoothing 
length, in units of the gravitational softening, /ism- The choice of 
this smoothing length is somewhat arbitrary. In radiative simula- 
tions, cooling causes gas to collapse and form high density cores 
at the centre of halos, thus causing the SPH smoothing length of 
cooled gas particles to fall below the gravitational softening. One 
therefore often restricts the SPH smoothing length to a certain frac- 
tion of the gravitational softening, ranging between zero (no restric- 
tion) and one (restricted to the gravitational softening). While the 
reference value adopted for our simulations is ftsm = 0.1, we have 
tested the effect of increasing it up to a value of unity. 

To separate the effects of numerics from those related to the 
physics included in the simulations, we decided to perform our 
study of numerical stability for the simplest version of the radiative 
runs, i.e. by not using the kinetic feedback associated with galactic 
winds {csfnw runs). 

In the left panel of Figure A2, we show the results of our reso- 
lution study. They confirm that the mass range over which the total 
subhalo mass function is reliably described widens from the Rl to 
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Figure A2. Comparison of the total (left panel) and stellar (right panel) subhalo mass functions, for runs with increasing resolution as well as for runs with 
the alternative 'gasification' scheme (see text for details). Whereas the total subhalo mass function seems to be converged for subhalos which are resolved by 
30-50 dark matter particles, the stellar subhalo mass function appears to converge only much more slowly. Here the one over eight gasification scheme seems 
to work better, because it converges earlier for fewer gas particles and only mildly more dark matter particles. 
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55 
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4.9e7 


5el0 
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Opteron 2.8 


R3 


1.5e8 
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6.6e8 


7.9e8 


4.0e8 


7el0 


63 
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1.5e8 


1.8e8 


8.8e7 


2el0 


257 


Opteron 2.8 



cooling in small subhalos, without overly increasing the CPU time 
requirements. 

Finally, Figure A3 shows the effect of changing the minimum 
value of the SPH smoothing length, /ism- As for the total subhalo 
mass function (right panel), it has only a mild effect at the smallest 
resolved masses. We note that no differences exist for hsm ^ 0.5, 
while a significant reduction of the mass function takes place for 
hs-m = 1- The effect is more apparent for the stellar mass func- 
tion (left panel). This result demonstrates that not allowing the SPH 
smoothing length to become smaller than the gravitational soften- 
ing makes the star clumps more fragile within the strong tidal field 
of the cluster. We note that using a small /ism is computationally 
cheaper, since it prevents the occurrence of computationally ex- 
pensive SPH computations for a large number of neighbours in the 
highly dense regions of cold gas. 



APPENDIX B: DATABASE 

The results of our substructure analysis for the dm, ovisc and 
csf version of all our simulated clusters (with the exception of 
g696) are publicly available within the German Virtual Observa- 
tory (GAVO). The access to the data is realized by a VO-oriented, 
SQL-queryable database similar to the one described in Lem- 
son & Virgo Consortium (2006). The web application that al- 
lows users to query the cluster database (called Hutt) is located 
at http : / / www . g-vo . or g/HydroClusters. Together with 
this publication, the data structures are made available to all users. 
The interface provides a full documentation of the available data 
and also gives different examples for scientific SQL-queries, return- 
ing data which can be used to produce similar plots than presented 
in this paper. In the future, we plan to allow users to register and 
to build up their own databases to store intermediate results from 
the queries to allow more complex analysis. Also, we plan to ex- 
tend the underlying database with further simulations that are not 
discussed in this study. 



Table Al. The masses of the different particle species for the runs at dif- 
ferent resolution. Column 1: run name. Columns 2-4: masses of the DM, 
gas and star particles, respectively (in units of Mq). Column 5: estimate of 
the stellar subhalo mass at which the simulations is numerically converged 
(since no higher resolution that R3 is achieved, we give only an upper limit 
in this case). Column 6: the CPU time (in hours) needed to perform the run. 



the R3 run. In general, runs at varying resolution provide similar 
results at masses corresponding to at least 20 DM particles, thus 
confirming the result shown in Fig. 3 for DM-only runs. Quite in- 
terestingly, the R2-8 and R3-8 runs, based on using equal masses 
for DM and gas particles, provide results virtually indistinguish- 
able from the corresponding R2 and R3 runs over the whole mass 
range where the results are stable against resolution. This confirms 
that the total mass function is determined by the DM mass resolu- 
tion and that a numerically converged total mass function can be 
obtained with the alternative way of "gasifying" initial conditions 
at a lower computational cost (see Column 6 in Table Al). 

The results are different for the stellar mass function, shown in 
the right panel of Fig. A2. Since these simulations are not includ- 
ing an effective feedback scheme able to regulate star formation, it 
is not surprising that the subhalo mass function progressively in- 
creases without any sign of convergence with resolution. A more 
interesting result concerns instead the behaviour of the R2-8 and 
R3-8 runs. Despite the fact that these runs have a poorer gas mass 
resolution than the RI and R2 runs, respectively, they produce a 
higher stellar subhalo mass function. For instance, the R3-8 run is 
consistent with the R3 run down to 1.5 x 10^" A/©, whereas the 
R2 run converges on the R3 run only at 5 x 10^" M0. This result 
is achieved with an acceptable increase of the computational cost. 
This highlights that increasing the DM mass resolution is indeed 
an efficient way of increasing the accuracy in the description of 




Figure A3. The total (left panel) and stellar (right panel) subhalo mass functions at different resolutions and using, at each resolution, different values for the 
minimum of the SPH smoothing length, hsm- 



